#include "gtest/gtest.h"
#include "sgSuperLU.h"

//#include "DataStructure/Common/include/Constant.h"

using namespace SG::Algebra;  // Open the namespace
// using SG::DataStructure::Common::TOLERANCE;

class SuperLUTest : public ::testing::Test
{
public:
    static void SetUpTestSuite ()
    {
        s_solver.reset(new sgSuperLU());
    }

    static std::shared_ptr<LinearSolver> s_solver;
};

std::shared_ptr<LinearSolver> SuperLUTest::s_solver;

/*
* Test the A*x = b given by SuperLU User Guide.
* where A is a 5X5 matrix and b is a 5X1 matrix.
* 
* A =
*   [
*       19.00,  00.00,  21.00,  21.00, 00.00;
*       12.00,  21.00,  00.00,  00.00, 00.00;
*       00.00,  12.00,  16.00,  00.00, 00.00;
*       00.00,  00.00,  00.00,  05.00, 21.00;
*       12.00,  12.00,  00.00,  00.00, 18.00  
*   ]
*
* b = [1.0; 1.0; 1.0; 1.0; 1.0]
*
* OctaveOnline give the solution as followings,
* x = inv(A)*b = [-0.031250; 0.065476; 0.013393; 0.062500; 0.032738]
*
* see https://octave-online.net/
*/
TEST_F (SuperLUTest, decomposTest)
{
    typedef Eigen::Triplet<double> Triplet;
    std::vector<Triplet> triplets;

    triplets.push_back(Triplet(0, 0, 19.0));
    triplets.push_back(Triplet(0, 2, 21.0));
    triplets.push_back(Triplet(0, 3, 21.0));

    triplets.push_back(Triplet(1, 0, 12.0));
    triplets.push_back(Triplet(1, 1, 21.0));

    triplets.push_back(Triplet(2, 1, 12.0));
    triplets.push_back(Triplet(2, 2, 16.0));

    triplets.push_back(Triplet(3, 3, 5.0));
    triplets.push_back(Triplet(3, 4, 21.0));

    triplets.push_back(Triplet(4, 0, 12.0));
    triplets.push_back(Triplet(4, 1, 12.0));
    triplets.push_back(Triplet(4, 4, 18.0));

    int m = 5; // row count
    int n = 5; // column count

    SpMat_t mat(m, n);
    mat.setFromTriplets(triplets.begin(), triplets.end());    

    //std::cout << std::endl; 
    //std::cout << "OuterStarts: ";
    //for (int i = 0; i < mat.outerSize() + 1; ++ i)
    //{
    //    std::cout << mat.outerIndexPtr()[i] << " ";
    //}
    //std::cout << std::endl ;

    //std::cout << "InnerIndices: ";
    // for (int i = 0; i < mat.nonZeros(); ++ i)
    // {
    //     //std::cout << mat.innerIndexPtr()[i] << " ";
    //     std::cout << mat.innerIndexPtr()[i] << " ";
    // }
    // std::cout  << std::endl;

    // std::cout << "Values: ";
    // for (int i = 0; i < mat.nonZeros(); ++ i)
    // {
    //     std::cout << mat.valuePtr()[i] << " ";
    // }
    // std::cout << std::endl;   
    // std::cout << std::endl;   

    s_solver->compute(mat);
    auto info = s_solver->info();    

    //Check
    ASSERT_TRUE (Info_t::Success == info);

}

TEST_F (SuperLUTest, solveTest)
{
    int m = 5; // row count;
    //int rhs = 1;

    // solution from OctaveOnline
    // x = inv(A)*b
    // see https://octave-online.net/
    double exact[] = {-0.031250, 0.065476, 0.013393, 0.062500, 0.032738};

    Vec_t b(m);
    for (int i = 0; i < m; ++ i )
    {
        b[i] = 1.0;
    }

    auto x = s_solver->solve(b);
    auto info = s_solver->info();

    //Check
    ASSERT_TRUE (Info_t::Success == info);   

    //for (int i = 0; i < m; ++ i)
    //{
    //    std::cout << exact[i] << ", " << x[i] << std::endl;
    //}

    for (int i = 0; i < m; ++ i)
    {
        ASSERT_NEAR (exact[i], x(i), TOLERANCE);
    }   
}